﻿namespace RayDen.RayEngine.IntersectionEngines.OpenCL {
    internal partial class Kernels {


        public static string LuxBVHKernel = @"
typedef struct {
	float x, y, z;
} Point;

typedef struct {
	float x, y, z;
} Vector;

typedef struct {
	Point o;
	Vector d;
	float mint, maxt;
} Ray;

typedef struct {
	float Distance;
	float U, V; // Barycentric coordinates of the hit point
	unsigned int Index;
} RayHit;

typedef struct {
	unsigned int v[3];
} Triangle;

typedef struct {
	Point pMin, pMax;
} BBox;

typedef struct {
	BBox bbox;
	unsigned int primitive;
	unsigned int skipIndex;
} GpuBVHNode;

void TriangleIntersect(
		const float4 rayOrig,
		const float4 rayDir,
		const float minT,
		float *maxT,
		unsigned int *hitIndex,
		float *hitB1,
		float *hitB2,
		const unsigned int currentIndex,
		__global Point *verts,
		__global Triangle *tris) {

	// Load triangle vertices
	__global Point *p0 = &verts[tris[currentIndex].v[0]];
	__global Point *p1 = &verts[tris[currentIndex].v[1]];
	__global Point *p2 = &verts[tris[currentIndex].v[2]];

	float4 v0 = (float4) (p0->x, p0->y, p0->z, 0.f);
	float4 v1 = (float4) (p1->x, p1->y, p1->z, 0.f);
	float4 v2 = (float4) (p2->x, p2->y, p2->z, 0.f);

	// Calculate intersection
	float4 e1 = v1 - v0;
	float4 e2 = v2 - v0;
	float4 s1 = cross(rayDir, e2);

	const float divisor = dot(s1, e1);
	if (divisor == 0.f)
		return;

	const float invDivisor = 1.f / divisor;

	// Compute first barycentric coordinate
	const float4 d = rayOrig - v0;
	const float b1 = dot(d, s1) * invDivisor;
	if (b1 < 0.f)
		return;

	// Compute second barycentric coordinate
	const float4 s2 = cross(d, e1);
	const float b2 = dot(rayDir, s2) * invDivisor;
	if (b2 < 0.f)
		return;

	const float b0 = 1.f - b1 - b2;
	if (b0 < 0.f)
		return;

	// Compute _t_ to intersection point
	const float t = dot(e2, s2) * invDivisor;
	if (t < minT || t > *maxT)
		return;

	*maxT = t;
	*hitB1 = b1;
	*hitB2 = b2;
	*hitIndex = currentIndex;
}

int BBoxIntersectP(
		const float4 rayOrig, const float4 invRayDir,
		const float mint, const float maxt,
		const float4 pMin, const float4 pMax) {
	const float4 l1 = (pMin - rayOrig) * invRayDir;
	const float4 l2 = (pMax - rayOrig) * invRayDir;
	const float4 tNear = fmin(l1, l2);
	const float4 tFar = fmax(l1, l2);

	float t0 = max(max(max(tNear.x, tNear.y), max(tNear.x, tNear.z)), mint);
    float t1 = min(min(min(tFar.x, tFar.y), min(tFar.x, tFar.z)), maxt);

	return (t1 > t0);
}

__kernel void Intersect(
		__global Ray *rays,
		__global RayHit *rayHits,
		__global Point *verts,
		__global Triangle *tris,
		const unsigned int triangleCount,
		const unsigned int nodeCount,
		__global GpuBVHNode *bvhTree,
		const unsigned int rayCount) {
	// Select the ray to check
	const int gid = get_global_id(0);
	if (gid >= rayCount)
		return;

	float4 rayOrig,rayDir;
	float minT, maxT;
	{
		__global float4 *basePtr =(__global float4 *)&rays[gid];
		float4 data0 = (*basePtr++);
		float4 data1 = (*basePtr);

		rayOrig = (float4)(data0.x, data0.y, data0.z, 0.f);
		rayDir = (float4)(data0.w, data1.x, data1.y, 0.f);

		minT = data1.z;
		maxT = data1.w;
	}

	//float4 rayOrig = (float4) (rays[gid].o.x, rays[gid].o.y, rays[gid].o.z, 0.f);
	//float4 rayDir = (float4) (rays[gid].d.x, rays[gid].d.y, rays[gid].d.z, 0.f);
	//float minT = rays[gid].mint;
	//float maxT = rays[gid].maxt;

	float4 invRayDir = (float4) 1.f / rayDir;

	unsigned int hitIndex = 0xffffffffu;
	unsigned int currentNode = 0; // Root Node
	float b1, b2;
	unsigned int stopNode = bvhTree[0].skipIndex; // Non-existent

	float4 pMin, pMax, data0, data1;
	__global float4 *basePtr;
	while (currentNode < stopNode) {

		basePtr =(__global float4 *)&bvhTree[currentNode];
		data0 = (*basePtr++);
		data1 = (*basePtr);

		pMin = (float4)(data0.x, data0.y, data0.z, 0.f);
		pMax = (float4)(data0.w, data1.x, data1.y, 0.f);

		if (BBoxIntersectP(rayOrig, invRayDir, minT, maxT, pMin, pMax)) {
			//const unsigned int triIndex = bvhTree[currentNode].primitive;
			const unsigned int triIndex = as_uint(data1.z);

			if (triIndex != 0xffffffffu)
				TriangleIntersect(rayOrig, rayDir, minT, &maxT, &hitIndex, &b1, &b2, triIndex, verts, tris);

			currentNode++;
		} else {
			//bvhTree[currentNode].skipIndex;
			currentNode = as_uint(data1.w);
		}
	}

	// Write result
	rayHits[gid].Distance   = maxT;
	rayHits[gid].U = b1;
	rayHits[gid].V = b2;
	rayHits[gid].Index = hitIndex;
} 
";



        #region SOA

        public static string SoaVHKernel = @"
typedef struct {
	float x, y, z;
} Point;

typedef struct {
	float x, y, z;
} Vector;

typedef struct {
	Point o;
	Vector d;
	float mint, maxt;
} Ray;

typedef struct {
	float Distance;
	float U, V; // Barycentric coordinates of the hit point
	unsigned int Index;
} RayHit;

typedef struct {
	unsigned int v[3];
} Triangle;

typedef struct {
	Point pMin, pMax;
} BBox;

typedef struct {
	BBox bbox;
	unsigned int primitive;
	unsigned int skipIndex;
} GpuBVHNode;

void TriangleIntersect(
		const float4 rayOrig,
		const float4 rayDir,
		const float minT,
		float *maxT,
		unsigned int *hitIndex,
		float *hitB1,
		float *hitB2,
		const unsigned int currentIndex,
		__global float *V0,
		__global float *V1,
		__global float *V2,
		__global uint *tv0,
		__global uint *tv1,
		__global uint *tv2) {

	// Load triangle vertices
	float3 p0 = (float3)(V0[tv0[currentIndex]],	V1[tv0[currentIndex]], V2[tv0[currentIndex]]);
	float3 p1 = (float3)(V0[tv1[currentIndex]],	V1[tv1[currentIndex]], V2[tv1[currentIndex]]);
	float3 p2 = (float3)(V0[tv2[currentIndex]],	V1[tv2[currentIndex]], V2[tv2[currentIndex]]);


	float4 v0 = (float4) (p0.x, p0.y, p0.z, 0.f);
	float4 v1 = (float4) (p1.x, p1.y, p1.z, 0.f);
	float4 v2 = (float4) (p2.x, p2.y, p2.z, 0.f);

	// Calculate intersection
	float4 e1 = v1 - v0;
	float4 e2 = v2 - v0;
	float4 s1 = cross(rayDir, e2);

	const float divisor = dot(s1, e1);
	if (divisor == 0.f)
		return;

	const float invDivisor = 1.f / divisor;

	// Compute first barycentric coordinate
	const float4 d = rayOrig - v0;
	const float b1 = dot(d, s1) * invDivisor;
	if (b1 < 0.f)
		return;

	// Compute second barycentric coordinate
	const float4 s2 = cross(d, e1);
	const float b2 = dot(rayDir, s2) * invDivisor;
	if (b2 < 0.f)
		return;

	const float b0 = 1.f - b1 - b2;
	if (b0 < 0.f)
		return;

	// Compute _t_ to intersection point
	const float t = dot(e2, s2) * invDivisor;
	if (t < minT || t > *maxT)
		return;

	*maxT = t;
	*hitB1 = b1;
	*hitB2 = b2;
	*hitIndex = currentIndex;
}

int BBoxIntersectP(
		const float4 rayOrig, const float4 invRayDir,
		const float mint, const float maxt,
		const float4 pMin, const float4 pMax) {
	const float4 l1 = (pMin - rayOrig) * invRayDir;
	const float4 l2 = (pMax - rayOrig) * invRayDir;
	const float4 tNear = fmin(l1, l2);
	const float4 tFar = fmax(l1, l2);

	float t0 = max(max(max(tNear.x, tNear.y), max(tNear.x, tNear.z)), mint);
    float t1 = min(min(min(tFar.x, tFar.y), min(tFar.x, tFar.z)), maxt);

	return (t1 > t0);
}
float4 ReadVecFromImg(int ind, __read_only image2d_t img)
{

   const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates
         CLK_ADDRESS_CLAMP | //Clamp to zeros
         CLK_FILTER_NEAREST; //Don't interpolate

   if (ind < 0) return 0;
   //Divide desired position by 4 because there are 4 components per pixel
   int imgPos = ind;

   int2 coords;
   coords.x = imgPos >> 13;
   coords.y = imgPos & 0x1fff;   //Reads the float4
   float4 temp = read_imagef(img, smp, coords);
    return temp;
}

uint4 ReadFromImg(int ind, __read_only image2d_t img)
{

   const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates
         CLK_ADDRESS_CLAMP | //Clamp to zeros
         CLK_FILTER_NEAREST; //Don't interpolate

   if (ind < 0) return 0;
   //Divide desired position by 4 because there are 4 components per pixel
   int imgPos = ind;

   int2 coords;
   coords.x = imgPos >> 13;
   coords.y = imgPos & 0x1fff;   //Reads the float4
   uint4 temp = read_imageui(img, smp, coords);
   return temp;
}
__kernel void Intersect(
		__global Ray *rays,
		__global RayHit *rayHits,
		__global float *v0,
		__global float *v1,
		__global float *v2,
		__global uint *tv0,
		__global uint *tv1,
		__global uint *tv2,
		const unsigned int triangleCount,
		const unsigned int nodeCount,
        __read_only  image2d_t bvhTree,
		//__global GpuBVHNode *bvhTree,
		const unsigned int rayCount) {
	// Select the ray to check
	const int gid = get_global_id(0);
	if (gid >= rayCount)
		return;

	float4 rayOrig,rayDir;
	float minT, maxT;
	{
		__global float4 *basePtr =(__global float4 *)&rays[gid];
		float4 data0 = (*basePtr++);
		float4 data1 = (*basePtr);

		rayOrig = (float4)(data0.x, data0.y, data0.z, 0.f);
		rayDir = (float4)(data0.w, data1.x, data1.y, 0.f);

		minT = data1.z;
		maxT = data1.w;
	}
    const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates
                      CLK_ADDRESS_CLAMP | //Clamp to zeros
                      CLK_FILTER_NEAREST;

	float4 invRayDir = (float4) 1.f / rayDir;

	unsigned int hitIndex = 0xffffffffu;
	unsigned int currentNode = 0; // Root Node
	float b1, b2;
    
	unsigned int stopNode = ReadFromImg(1, bvhTree).w; //bvhTree[0].skipIndex; // Non-existent

	float4 pMin, pMax, data0, data1;
	//__global float4 *basePtr;
	while (currentNode < stopNode) {

		//basePtr = (__global float4 *)&bvhTree[currentNode];
		data0 = ReadVecFromImg(currentNode, bvhTree);
		data1 = ReadVecFromImg(currentNode+1, bvhTree);

		pMin = (float4)(data0.x, data0.y, data0.z, 0.f);
		pMax = (float4)(data0.w, data1.x, data1.y, 0.f);

		if (BBoxIntersectP(rayOrig, invRayDir, minT, maxT, pMin, pMax)) {
			//const unsigned int triIndex = bvhTree[currentNode].primitive;
			const unsigned int triIndex = as_uint(data1.z);

			if (triIndex != 0xffffffffu)
				TriangleIntersect(rayOrig, rayDir, minT, &maxT, &hitIndex, &b1, &b2, triIndex, v0, v1, v2 , tv0, tv1, tv2);

			currentNode++;
		} else {
			//bvhTree[currentNode].skipIndex;
			currentNode =  ReadFromImg(as_uint(data1.w), bvhTree).w;
		}
	}

	// Write result
	rayHits[gid].Distance   = maxT;
	rayHits[gid].U = b1;
	rayHits[gid].V = b2;
	rayHits[gid].Index = hitIndex;
} 
";
        #endregion


        #region OTHER
        public static string SimpleBVHKernel = @"
typedef struct {
	float x, y, z;
} Point;

typedef struct {
	float x, y, z;
} Vector;

typedef struct {
	Point o;
	Vector d;
	float mint, maxt;
} Ray;

typedef struct {
	float t;
	float b1, b2; // Barycentric coordinates of the hit point
	unsigned int index;
} RayHit;

typedef struct {
	unsigned int v[3];
} Triangle;

typedef struct {
	Point pMin, pMax;
} BBox;


typedef struct {
	BBox bbox;
    unsigned int index;
    unsigned char count;
    unsigned char axis;
 } GpuBVHNode;

static int TriangleIntersect(
		const float4 rayOrig,
		const float4 rayDir,
		const float minT,
		float *maxT,
		unsigned int *hitIndex,
		float *hitB1,
		float *hitB2,
		const unsigned int currentIndex,
		__global Point *verts,
		__global Triangle *tris) {

	// Load triangle vertices
	__global Point *p0 = &verts[tris[currentIndex].v[0]];
	__global Point *p1 = &verts[tris[currentIndex].v[1]];
	__global Point *p2 = &verts[tris[currentIndex].v[2]];

	float4 v0 = (float4) (p0->x, p0->y, p0->z, 0.f);
	float4 v1 = (float4) (p1->x, p1->y, p1->z, 0.f);
	float4 v2 = (float4) (p2->x, p2->y, p2->z, 0.f);

	// Calculate intersection
	float4 e1 = v1 - v0;
	float4 e2 = v2 - v0;
	float4 s1 = cross(rayDir, e2);

	const float divisor = dot(s1, e1);
	if (divisor == 0.f)
		return 0;

	const float invDivisor = 1.f / divisor;

	// Compute first barycentric coordinate
	const float4 d = rayOrig - v0;
	const float b1 = dot(d, s1) * invDivisor;
	if (b1 < 0.f)
		return 0;

	// Compute second barycentric coordinate
	const float4 s2 = cross(d, e1);
	const float b2 = dot(rayDir, s2) * invDivisor;
	if (b2 < 0.f)
		return 0;

	const float b0 = 1.f - b1 - b2;
	if (b0 < 0.f)
		return 0;

	// Compute _t_ to intersection point
	const float t = dot(e2, s2) * invDivisor;
	if (t < minT || t > *maxT)
		return 0;

	*maxT = t;
	*hitB1 = b1;
	*hitB2 = b2;
	*hitIndex = currentIndex;
    return 1;
}

static int BBoxIntersectP(
		const float4 rayOrig, const float4 invRayDir,
		const float mint, const float maxt,
		const float4 pMin, const float4 pMax) {
	const float4 l1 = (pMin - rayOrig) * invRayDir;
	const float4 l2 = (pMax - rayOrig) * invRayDir;
	const float4 tNear = fmin(l1, l2);
	const float4 tFar = fmax(l1, l2);

	float t0 = max(max(max(tNear.x, tNear.y), max(tNear.x, tNear.z)), mint);
    float t1 = min(min(min(tFar.x, tFar.y), min(tFar.x, tFar.z)), maxt);

	return (t1 > t0);
}

__kernel void Intersect(
		__global Ray *rays,
		__global RayHit *rayHits,
		__global Point *verts,
		__global Triangle *tris,
		const unsigned int triangleCount,
		const unsigned int nodeCount,
		__global GpuBVHNode *bvhTree,
		const unsigned int rayCount) {
	// Select the ray to check
	const int gid = get_global_id(0);
	if (gid >= rayCount)
		return;


	float4 rayOrig = (float4) (rays[gid].o.x, rays[gid].o.y, rays[gid].o.z, 0.f);
	float4 rayDir = (float4) (rays[gid].d.x, rays[gid].d.y, rays[gid].d.z, 0.f);
	float minT = rays[gid].mint;
	float maxT = rays[gid].maxt;

	float4 invRayDir = (float4) 1.f / rayDir;

	unsigned int hitIndex = 0xffffffffu;
	unsigned int currentNode = 0; // Root Node
	float b1 , b2, dist = 1e6f;
	unsigned int stopNode = bvhTree[0].index; // Non-existent

	float4 pMin, pMax, data0, data1;
    int todoOffset = 0;
    uint todo[32];
    int isDirNeg[3];

    isDirNeg[0] = invRayDir.x < 0.f;
    isDirNeg[1] = invRayDir.y < 0.f;
    isDirNeg[2] = invRayDir.z < 0.f;

	while (true) {
		float4 pMin = (float4)(bvhTree[currentNode].bbox.pMin.x,
				bvhTree[currentNode].bbox.pMin.y,
				bvhTree[currentNode].bbox.pMin.z,
				0.f);
		float4 pMax = (float4)(bvhTree[currentNode].bbox.pMax.x,
				bvhTree[currentNode].bbox.pMax.y,
				bvhTree[currentNode].bbox.pMax.z,
				0.f);
		if (BBoxIntersectP(rayOrig, invRayDir, minT, maxT, pMin, pMax)) {
			const unsigned int triIndex = bvhTree[currentNode].index;
			if (bvhTree[currentNode].count > 0 && TriangleIntersect(rayOrig, rayDir, minT, &maxT, &hitIndex, &b1, &b2, triIndex, verts, tris) && dist > maxT) {
                dist = maxT;
            } 
            else {
                if (isDirNeg[bvhTree[currentNode].axis] < 0.f) 
                {
                    todo[todoOffset++] = currentNode + 1;
                    currentNode = triIndex;
                }
                else 
                {
                    todo[todoOffset++] = triIndex;
                    currentNode = currentNode + 1;
                }
            }
            if (todoOffset == 0)
                break;
            currentNode = todo[--todoOffset];
		} 
        else 
        {
            if (todoOffset == 0)
                break;
            currentNode = todo[--todoOffset];
        }
	}

	// Write result
	rayHits[gid].t = maxT;
	rayHits[gid].b1 = b1;
	rayHits[gid].b2 = b2;
	rayHits[gid].index = hitIndex;
} 


__kernel void BFIntersect(
		__global Ray *rays,
		__global RayHit *rayHits,
		__global Point *verts,
		__global Triangle *tris,
		const unsigned int triangleCount) {
	// Select the ray to check
	const int gid = get_global_id(0);

	float4 rayOrig,rayDir;
	float minT, maxT;
	{
		__global float4 *basePtr =(__global float4 *)&rays[gid];
		float4 data0 = (*basePtr++);
		float4 data1 = (*basePtr);

		rayOrig = (float4)(data0.x, data0.y, data0.z, 0.f);
		rayDir = (float4)(data0.w, data1.x, data1.y, 0.f);

		minT = data1.z;
		maxT = data1.w;
	}

	//float4 rayOrig = (float4) (rays[gid].o.x, rays[gid].o.y, rays[gid].o.z, 0.f);
	//float4 rayDir = (float4) (rays[gid].d.x, rays[gid].d.y, rays[gid].d.z, 0.f);
	//float minT = rays[gid].mint;
	//float maxT = rays[gid].maxt;

	unsigned int hitIndex = 0xffffffffu, h1 = 0u;
	unsigned int currentNode = 0; // Root Node
	float b1, b2, dist = maxT, u,v;
	unsigned int stopNode = triangleCount;

	while (currentNode < stopNode) {
		const unsigned int triIndex = as_uint(currentNode);
		TriangleIntersect(rayOrig, rayDir, minT, &maxT, &hitIndex, &b1, &b2, triIndex, verts, tris);

        if (dist > maxT)
        {
            dist = maxT;
            u = b1;
            v = b2;
        }
        h1 = hitIndex;
		currentNode++;
	}

	rayHits[gid].t = dist;
	rayHits[gid].b1 = u;
	rayHits[gid].b2 = v;
	rayHits[gid].index = h1;
} 
        ";
        #endregion
    }
}
